// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

#ifndef SIMPLE_XYZ_H
#define SIMPLE_XYZ_H

// standard headers
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <vector>
#include <list>
#include <set>
#include <map>
#include <iomanip>
#include <cmath>
#include <cstdarg>
#include <stdexcept>
#include <utility>
#include <tuple>

#include <Eigen/Eigen>

/*
void load(std::istream & is, size_t & n, std::vector<Eigen::Vector3d> & lattice, std::vector<std::string> & labels, std::vector<Eigen::Vector3d> & coords, std::string & cline) 
{
	std::string line;
	if(std::getline(is,line))
	{
	        
		std::istringstream iss(line);
		iss >> n;
	}
        labels.resize(n);
        coords.resize(n);
        lattice.resize(3);
	if(std::getline(is,line))
	{
                cline = line;
		std::istringstream iss(line);
		iss >> lattice[0].x() >> lattice[0].y() >> lattice[0].z();
		iss >> lattice[1].x() >> lattice[1].y() >> lattice[1].z();
		iss >> lattice[2].x() >> lattice[2].y() >> lattice[2].z();
                // error check and bounding box
	}
	for(size_t i = 0; i < n; i++)
	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		iss >> labels[i];
		iss >> coords[i].x() >> coords[i].y() >> coords[i].z();
	}
        else throw std::runtime_error("Bad xyz file.");
};
*/

bool load(std::istream & is, size_t & n, std::vector<Eigen::Vector3d> & lattice, std::vector<std::string> & labels, std::vector<Eigen::Vector3d> & coords, std::string & cline) 
{
	std::string line;
	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		if(!(iss >> n)) return false;
	}
	else return false;
	
        labels.resize(n);
        coords.resize(n);
        lattice.resize(3);
	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		iss >> lattice[0].x() >> lattice[0].y() >> lattice[0].z();
		iss >> lattice[1].x() >> lattice[1].y() >> lattice[1].z();
		iss >> lattice[2].x() >> lattice[2].y() >> lattice[2].z();
		iss >> cline;
		// error check and bounding box
	}
	for(size_t i = 0; i < n; i++)
	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		iss >> labels[i];
		iss >> coords[i].x() >> coords[i].y() >> coords[i].z();
	}
        else throw std::runtime_error("Bad xyz file.");
        return true;
};


bool load(std::istream & is, size_t & n, std::vector<Eigen::Vector3d> & lattice, std::vector<std::string> & labels, std::vector<Eigen::Vector3d> & coords) 
{
	std::string line;
	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		if(!(iss >> n)) return false;
	}
	else return false;
	
        labels.resize(n);
        coords.resize(n);
        lattice.resize(3);
	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		iss >> lattice[0].x() >> lattice[0].y() >> lattice[0].z();
		iss >> lattice[1].x() >> lattice[1].y() >> lattice[1].z();
		iss >> lattice[2].x() >> lattice[2].y() >> lattice[2].z();
                // error check and bounding box
	}
	for(size_t i = 0; i < n; i++)
	if(std::getline(is,line))
	{
		std::istringstream iss(line);
		iss >> labels[i];
		iss >> coords[i].x() >> coords[i].y() >> coords[i].z();
	}
        else throw std::runtime_error("Bad xyz file.");
        return true;
};

void save(std::ostream & os,const size_t & n, const std::vector<Eigen::Vector3d> & lattice, const std::vector<std::string> & labels, const std::vector<Eigen::Vector3d> & coords) 
{
	os << std::left << std::setw(5) << n;
	os << std::endl;

	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << lattice[0].x() << ' ';
	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << lattice[0].y() << ' ';
	os << std::left << std::setw(16) << std::setprecision(6) << std::fixed << lattice[0].z() << ' ';

	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << lattice[1].x() << ' ';
	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << lattice[1].y() << ' ';
	os << std::left << std::setw(16) << std::setprecision(6) << std::fixed << lattice[1].z() << ' ';

	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << lattice[2].x() << ' ';
	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << lattice[2].y() << ' ';
	os << std::left << std::setw(12) << std::setprecision(6) << std::fixed << lattice[2].z() << ' ';
        
        os << std::endl;

	for(size_t i = 0; i < n; i++)
	{
		os << std::left << std::setw(5) << labels[i];

		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << coords[i].x() << ' ';
		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << coords[i].y() << ' ';
		os << std::right << std::setw(12) << std::setprecision(6) << std::fixed << coords[i].z() << ' ';
		os << std::endl;
	}
}

void print(std::ostream & os,const size_t & n, const std::vector<Eigen::Vector3d> & lattice, const std::vector<std::string> & labels, const std::vector<Eigen::Vector3d> & coords, const std::vector<bool> & selected) 
{
	os << std::left << std::setw(5) << n;
	os << std::endl;

	os << std::left << std::setw(6) << std::setprecision(3) << std::fixed << lattice[0].x() << ' ';
	os << std::left << std::setw(6) << std::setprecision(3) << std::fixed << lattice[0].y() << ' ';
	os << std::left << std::setw(8) << std::setprecision(3) << std::fixed << lattice[0].z() << ' ';

	os << std::left << std::setw(6) << std::setprecision(3) << std::fixed << lattice[1].x() << ' ';
	os << std::left << std::setw(6) << std::setprecision(3) << std::fixed << lattice[1].y() << ' ';
	os << std::left << std::setw(8) << std::setprecision(3) << std::fixed << lattice[1].z() << ' ';

	os << std::left << std::setw(6) << std::setprecision(3) << std::fixed << lattice[2].x() << ' ';
	os << std::left << std::setw(6) << std::setprecision(3) << std::fixed << lattice[2].y() << ' ';
	os << std::left << std::setw(8) << std::setprecision(3) << std::fixed << lattice[2].z() << ' ';
        
        os << std::endl;

	for(size_t i = 0; i < n; i++)
	{
		os << std::left << std::setw(5) << labels[i];

		os << std::right << std::setw(6) << std::setprecision(3) << std::fixed << coords[i].x() << ' ';
		os << std::right << std::setw(6) << std::setprecision(3) << std::fixed << coords[i].y() << ' ';
		os << std::right << std::setw(6) << std::setprecision(3) << std::fixed << coords[i].z() << ' ';
                os << std::right << std::setw(6) << std::setprecision(3) << std::fixed << i+1 << ' ';
                os << std::right << std::setw(6) << std::setprecision(3) << std::fixed << int(i)-int(n) << ' ';
                
		os << std::endl;
	}
}


inline int convert_index(const size_t & n,int i)
{
        i =  i < 0 ? int(n) + i : i - 1; 
        if( (i<0) || (i >= int(n)) ) throw std::runtime_error("Wrong atom index!");
        return  i; 
}

struct xyz
{
        size_t n;
        std::vector<Eigen::Vector3d> lattice;
        std::vector<std::string> labels;
        std::vector<Eigen::Vector3d> coords;
        std::vector<bool> selected;
        std::string comment;
};

void load(std::istream & is,xyz & s)
{
        load(is,s.n,s.lattice,s.labels,s.coords,s.comment);
        s.selected.resize(s.n); //all false?
}

void load(std::string fn,xyz & s)
{
	std::ifstream is(fn.c_str());
        load(is,s);
	is.close();
}

void load(std::string fn,std::vector<xyz> & sv)
{
        sv.clear();
	std::ifstream is(fn.c_str());
        xyz s;
        while(load(is,s.n,s.lattice,s.labels,s.coords,s.comment))
        {
            s.selected.resize(s.n); //all false?
            sv.push_back(s);
        }
	is.close();
}

void save(std::ostream & os,const xyz & s)
{
        save(os,s.n,s.lattice,s.labels,s.coords);
}
void save(std::string fn,const xyz & s)
{
        std::ofstream os(fn.c_str());
        save(os,s);
	os.close();
}

void save(std::string fn,const std::vector<xyz> & sv)
{
	std::ofstream os(fn.c_str());
	for(size_t i = 0;i<sv.size();i++)
	save(os,sv[i]);
	os.close();
}


void print(std::ostream & os,const xyz & s)
{
        print(os,s.n,s.lattice,s.labels,s.coords,s.selected);
}

xyz & select(xyz & s,const std::list<int> & is)
{
        for(std::list<int>::const_iterator it = is.begin(); it != is.end(); it++)
        s.selected[convert_index(s.n,*it)] = true;
        return s;
}


xyz & neg(xyz & s)
{
        for(size_t i = 0; i < s.n; i++) s.selected[i] = !s.selected[i];
        return s;
}



xyz & operator += (xyz & s,Eigen::Vector3d & v)
{
        for(size_t i = 0; i < s.n; i++)
        if(s.selected[i])
        s.coords[i] += v;
        return s;
}

xyz & operator -= (xyz & s,Eigen::Vector3d & v)
{
        for(size_t i = 0; i < s.n; i++)
        if(s.selected[i])
        s.coords[i] -= v;
        return s;
}

/*
xyz & apply(xyz & s,Eigen::Matrix3d & m)
{
        for(size_t i = 0; i < s.n; i++)
        if(s.selected[i])
        s.coords[i] = m*s.coords[i];
        return s;
}
*/

xyz & apply(xyz & s,Eigen::Transform<double,3,Eigen::Affine> & m)
{
        for(size_t i = 0; i < s.n; i++)
        if(s.selected[i])
        s.coords[i] = m*s.coords[i];
        return s;
}

xyz & pour(xyz & s,const xyz & s2)
{
        for(size_t i = 0; i < s2.n; i++)
        if(s2.selected[i])
        {       
                s.n++;
                s.coords.push_back(s2.coords[i]);
                s.labels.push_back(s2.labels[i]);
                s.selected.push_back(true);
        }
        return s;
}

xyz & sort(xyz & s, bool (*greater)(const Eigen::Vector3d&,const Eigen::Vector3d&,double),double tiny )
{
        for(size_t i = 0; i < s.n; i++)
        if(s.selected[i])
        {
            for(size_t j = i+1; j < s.n; j++)
            if(s.selected[j])
            if(greater(s.coords[i],s.coords[j],tiny)) 
            {
                std::swap(s.coords[j],s.coords[i]);
                std::swap(s.labels[j],s.labels[i]);
            }
        }
      
}


/* how this is compatible with selected ?
xyz & reverse_order(xyz & s)
{
        for(size_t i = 0; i < s.n; i++)
        if(s.selected[i])
        {
            size_t minj = i;
            for(size_t j = i+1; j < s.n; j++)
            if(s.selected[j])
            if(greater(s.coords[i],s.coords[j])) 
            {
                minj = j;
                if(minj != i)
                std::swap(s.coords[minj],s.coords[i]);
            }
        }
      
}
*/

void bbox(xyz & s,Eigen::Vector3d & bboxmin,Eigen::Vector3d & bboxmax)
{
    bboxmin = s.coords[0];                                                                                                                                                        
    bboxmax = s.coords[0];                                                                                                                                                        
    for(int i=0; i<s.n; i++)                                                                                                                                                        
    if(s.selected[i])
    {                                                                                                                                                                             
        if(s.coords[i].x() < bboxmin.x()) bboxmin.x() = s.coords[i].x();                                                                                                              
        if(s.coords[i].y() < bboxmin.y()) bboxmin.y() = s.coords[i].y();                                                                                                              
        if(s.coords[i].z() < bboxmin.z()) bboxmin.z() = s.coords[i].z();                                                                                                              
        if(s.coords[i].x() > bboxmax.x()) bboxmax.x() = s.coords[i].x();                                                                                                              
        if(s.coords[i].y() > bboxmax.y()) bboxmax.y() = s.coords[i].y();                                                                                                              
        if(s.coords[i].z() > bboxmax.z()) bboxmax.z() = s.coords[i].z();                                                                                                              
    }                                                                                                                                                                             
}                                                                                                                                                                                                                                                                                                                                                                                                  

                                                                                                                                                                                                                                                                                                              
#endif // SIMPLE_XYZ_H


